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ABSTRACT 



An Implicit Fi ni te -Dif f erence ( IFD) program that 
incorporates exact interface conditions has been developed 
for solving the parabolic equation. The model preserves 
continuity of pressure and continuity of the normal 
component of particle velocity at the interface between 
media having different sound speeds and densities. 
Interface conditions are preserved for horizontal and 
sloping interfaces along a user-specified bottom profile. 
Test cases are included to demonstrate the use of the model. 
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I. 



INTRODUCTION 



Since its introduction to the underwater acoustics 
community (Hardin and Tappert, 1973), the parabolic wave 
equation has stimulated a considerable amount of interest. 
The first solution programs used a split-step fast Fourier 
transform method to solve the parabolic equation; however, 
other solution techniques have been developed (McDaniel and 
Lee, 1982). One of the motives for developing alternative 
solution techniques is that problems arise when the Fourier 
transform encounters an interface between two media having 
different sound speeds or densities (Lee and Botseas, 1982). 

One alternative solution technique uses an implicit 
finite-difference (IFD) solution method. The IFD method is 
unconditionally stable and has the capability to incorporate 
desired interface conditions. Implicit finite-difference 
methods for solving parabolic equations have been studied 
extensively by many authors. 

A computer program that utilizes the IFD method to solve 
the parabolic equation has been developed and is examined in 
detail in this thesis. The computer program predicts 
acoustic propagation loss in environments having user- 
specified bottom profiles. The program preserves continuity 
of pressure and continuity of the normal component of 
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particle velocity at an interface between media having 
different sound speeds and densities. 

The program utilizes concepts developed by earlier 
authors. The use of the IFD method to solve the parabolic 
equation in underwater acoustics was developed by Lee and 
Papadakis (1979). The mathematical treatment of horizontal 
and sloping interfaces was developed by McDaniel and Lee 
(1 982) and Lee and McDaniel (1 983). And finally, the 
program utilizes some design features of an earlier computer 
program developed by Lee and Botseas (1982). 



II. PARABOLIC EQUATION 



A. INTRODUCTION 

The parabolic equation is an approximation to the 
elliptical wave equation. The derivation of the parabolic 
equation begins with the reduced wave equation (Helmholtz 
equation) in the form 

V 2 P + k 2 = 0 2.1 

or 

V 2 P + k 0 2n2 P = 0 2.2 

where 

k = wave number ( = w/ c ) 

k Q = reference wave number (= w/c Q ) 

n = index of refraction (= c Q /c) 

p = time independent factor of complex pressure 

c = sound speed 

c Q = reference sound speed 

w = angular source frequency (= 2 7T'f) 

For the case of cylindrical symmetry (2.2) becomes 

p rr + (1/r) p r + p zz + k 0 2 n 2 p =0 2.3 

It is then assumed that p is of the form 

p = u(r , z ) S(r) 

where u is a function of both range and depth and S is a 
function of range only. Substitution of (2.4) into (2.3) 
and separation of variables shows that S(r) must satisfy 



Bessel's equation of zero-order. For exp (-iwt) time 
dependence and outgoing waves, the solution is the zeroth- 
order Hankel function of the first kind, 

S(r) = H q ( 1 ) ( k Q r ) . 

Further, u(r,z) must satisfy 



1 2 2 

u rr + u zz + + - S r ) u r + k Q (n^ - 1 ) u = 0. 

r S 



2.5 



With the help of the far-field asymptotic approximation for 
the Hankel function, (2.5) can be reduced to 



u rr + u zz 



+ 2ik rt u_ + k-/(n^ - 1) u = 0. 



L o“r ' ' v o 

We now assume that u varies slowly with respect to range, 



2.6 



u 



rr 



< < 



2 k o u r 



combining (2.6) and (2.7) results in 



u 



zz 



+ 2ik Q u r + k Q 2 (n 2 - 1 ) u = 0. 



2.7 



2.8 



Rearranging (2.8) results in the parabolic equation in the 
form 



u = a(k ,r,z) u + b(k ,r,z) u 



zz 



2.9 



where 

a(k Q ,r,z) = (ik Q /2) [n 2 (r,z) - 1] 

b(k Q ,r,z) = i/2k 0 . 

The assumption (2.7), fundamental to the parabolic equation 
method, is equivalent to neglecting back-scatter. 



B. SPLIT-STEP FAST FOURIER TRANSFORM SOLUTION 
1 . Description 

For the first few years after its introduction into 
the acoustical community the parabolic equation was solved 
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exclusively with the help of the split-step fast Fourier 
transform (SSFFT) method developed by Tappert and Hardin 
(Jensen and Krol, 1975). In this method, u 2z in (2.8) is 
represented by the inverse transform of its Fourier 
transform. The SSFFT method requires periodic boundary 
conditions in z because of the finite Fourier transform. 
This is handled by introducing an artificial, horizontal, 
pressure release bottom below the physical bottom. The 
SSFFT method is unconditionally stable (Brock, 1978). 

The SSFFT method has been implemented by Jensen and 
Krol and by Brock. Detailed descriptions can be found in 
publications of Jensen and Krol (1975) and Brock (1978). 

2 . Interface Treatment 

Errors introduced by the SSFFT method are 
proportional to the range step and to n zz where n is the 
index of refraction (Jensen and Krol, 1975). Because the 
index of refraction has a large change across the water- 
sediment interface n zz will be large and thus the error will 
be large. To reduce this error, a very small horizontal 
range step must be used. However, this results in very long 
computer execution time. The problem of a discontinuity in 
sound speed at the water-sediment interface and the 
resultant difficulties in solving shallow water propagation 
problems using the SSFFT method are addressed in Jensen and 
Krol (1975). 
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Another more serious problem with the SSFFT method 
is that it neglects any density difference between the water 
and the sediment. A density difference can be important in 
that it influences the reflection coefficient. The problem 
becomes more significant as the density difference becomes 
larger. 

In summary, the discontinuities in sound speed and 
in density at the wa ter- sediment interface cause problems 
for the SSFFT method. The SSFFT method is therefore 
intrinsically better suited for deep water propagation 
environments for which the water -sediment interface is an 
unimportant feature. 

C. IMPLICIT FINITE-DIFFERENCE SOLUTION METHOD 

In 1 979 Lee and Papadakis introduced the Crank-Nicolson 
implicit finite-difference method to solve the parabolic 
equation for underwater acoustic propagation. The Crank- 
Nicolson method uses a second-order central difference 
formula to approximate u 2Z in (2.9) and casts the problem in 
the form of a tridiagonal matrix system. A representative 
row in the matrix system (the m 1 -* 1 row) is 
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2.10 



where 

k = horizontal range increment 

h = vertical depth increment 

Supercripts are used to indicate range indices and 
subscripts are used to indicate depth indices. In (2.10) 
the field is known at range index n and is to be solved at 
range index n + 1. Therefore, the right hand side of (2.10) 
reduces to a single, known value and the solution field is 
advanced from range index n to range index n+1 by solving 
the tridiagonal system of equations. 

The IFD scheme is consistant, unconditionally stable and 
it converges to the theoretical solution as the range and 
depth increments tend to zero (Lee et al., 1981). An 
advantage of selecting an implicit scheme over an explicit 
scheme is that an explicit scheme is only conditionally 
stable (Lee and Papadakis, 1979). Another advantage of the 
implicit scheme is smaller errors. More detailed 
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information addressing both implicit and explicit solutions 
of parabolic equations can be found in Gerald (1980). 

The first IFD scheme handled discontinuities in the 
speed of sound profile but did not consider the effects of 
density discontinuities. It therefore did not correctly 
treat the interface between media having different 
densities . 

D. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 

HORIZONTAL INTERFACE 

In 1982 McDaniel and Lee introduced a scheme for 

\ 

incorporating a horizontal interface into the IFD method. 
The interface separates two media with different sound 
speeds and densities (Figure 1). 




Figure 1 . IFD Treatment of a Horizontal Interface 



16 



The scheme preserves continuity of pressure and continuity 
of the normal component of particle velocity across the 
interface and does not affect the stability of the IFD 
method. 

The interface equation in the tridiagonal matrix system 
that results from incorporating continuity of pressure and 
continuity of the normal component of particle velocity is 
(McDaniel and Lee, 1982) 



S P n+1 , l-U P n+1 Q n+1 
^2 m 2 m m 





h 



2.11 




k ~n 






/fl 

^2 



k p n 



u 



n 

L 

m 



<1 + ZT-), 



h 2 ”^2 




where 




P 



b 



a 1 a 2 



Q 



+ 



b 1 ^2 b 2 
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/"l = density in layer 1 (water) 

/°2 ~ density in layer 2 (sediment) 
a and b are defined in (2.9) 

Incorporating the horizontal interface into the IFD method 
requires inserting (2.10) for the row in the tridiagonal 
matrix system that corresponds to the interface. 

The error in the solution is 
0 (k 3 + kh) 
on the interface and 
O (k 3 + kh 2 ) 

in a continuous medium (McDaniel and Lee, 1982). 

E. IMPLICIT FINITE-DIFFERENCE METHOD: TREATMENT OF A 

SLOPING INTERFACE 

In 1983, Lee and McDaniel extended their treatment of an 
interface between two media to include the case of a sloping 
interface. As for the case of a horizontal interface, the 
treatment of a sloping interface preserves continuity of 
pressure and continuity of the normal component of particle 
velocity at the interface between media having different 
densities and sound speeds. 

The problem of a sloping interface is separated into two 
cases: downslope, and upslope. Each case requires inserting 
two new rows into the IFD tridiagonal matrix system. One 
new row is required at the level corresponding to the 
interface at the range where the solution is known and the 
second new row is required at the level corresponding to the 
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interface at the range at which the solution is to be 
solved. Therefore, four new equations are required to 
cover both the downslope and upslope cases. The four 
sloping interface equations are derived and shown in Lee and 
McDaniel (1983). The equations are somewhat involved but 
they are of the same tridiagonal form as the original IFD 
matrix equations. The error for points on or adjacent to a 
sloping interface is (Lee and McDaniel, 1983) 

0 (k^ + kh) . 



19 



III. COMPUTER IMPLEMENTATION 



A. INTRODUCTION 

An IFD solution program that implements (2.10), (2.11 ) 
and the four sloping interface equations has been developed 
to solve underwater propagation problems. The program is 
written in FORTRAN using single precision, complex 
arithmetic and has been installed on the Naval Postgraduate 
School's IBM-3033 digital computer. Appendix A contains a 
program listing. 

The solution program consists of one main program and 20 
subroutines. A modular program construction was selected 
for flexibility and clarity. 

Within each routine a structured programming approach is 
utilized. The structured program format, coupled with 
generous commenting, makes the program relatively easy to 
trace through. 

As presently installed on the NPS computer the solution 
program is run interactively. Appendix B contains details. 

B. GENERAL CHARACTERISTICS 

The solution program handles the following environmental 
conditions: range independent sound speed profile in the 
water column, range dependent bottom profile and iso-speed, 
iso-density sedimentary bottom layer. The program utilizes 
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a Gaussian starting field and an artificial pressure release 
surface in the sediment at a user-specified depth. (An 
artificial pressure release surface is not required for the 
IFD solution method; however, such a surface permits 
efficient solution for the pressure field when it is known 
to tend to zero at great depths in the bottom.) 

Attenuation in the water and sediment is introduced 
using complex indices of refraction. Artificially strong 
attenuation is applied in the lower portion of the sediment 
layer to enhance attenuation of the field above the 
artificial pressure release surface. 

A user-specified bottom profile is input as a series of 
linear bottom segments. The range step along a horizontal 
bottom segment is set to a user-specified value. The range 
step along a sloping bottom segment is automatically set by 
the program so that the sloping bottom intersects the next 
vertical grid point. Bottom modifications are required in 
certain situations to meet the requirements that the range 
step not be too large and that the interface must pass 
through a grid point at every range at which the pressure 
field is solved. For a very gently sloping bottom, if the 
calculated range step exceeds a user- spec if ied value then 
the program will automatically model the bottom as a series 
of level and sloping sections. For these cases, the 
difference between the modified bottom and the user- 
specified bottom is always less than or equal to one-half 
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the vertical grid increment. The user is informed if the 
bottom is modified. 

It is foreseen that future enhancements will increase 
the program's generality. In particular, changes to allow 
range dependent sound speed profiles, sound speed profiles 
in the sediment and a user-specified starting field should 
be relatively simple. The modular construction of the 
program facilitates these types of changes. 

C. MAIN PROGRAM IFD 

IFD is the main program. It controls program execution 
and calls subroutines as appropriate. 

The first executable statement in IFD calls subroutine 
ERRSET, a system subroutine peculiar to the NPS computer 
that correctly sets a variable value to zero when an 
underflow condition exists. Most computer systems do this 
automatically; however, depending on the particular system a 
call similar to ERRSET may be required. 

IFD then calls subroutine READ to read input data, 
subroutine SVPW to calculate the sound speed at grid points 
in the water column, subroutine INITAL to initialize 
constants and variables, and then subroutine MATCON to 
calculate matrix constants. Subroutine SFIELD is then 
called to calculate the Gaussian starting field, followed by 
subroutines WRITE1 and PRINT1 which write and print output 
data. In the context of this program, "write" refers to 
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writing unformatted data into a file to be used by the 
plotting routine and "print" refers to writing formatted 
data into a file which can be sent directly to the printer. 

IFD then calls subroutine NEWSEG which is the beginning 
of a loop that is called every time a new linear bottom 
segment is reached. NEWSEG calculates variables that 
characterize a new bottom segment. The next call is to 
subroutine NEWMAT which calculates matrix elements for the 
new bottom segment and advances the solution field one range 
step. IFD then enters a loop that advances the solution one 
range step for every pass through the loop. Inside the loop 
the range markers are advanced and the solution is advanced 
one step for the downslope, level, upslope, modified bottom 
downslope or modified bottom upslope situation as 
appropriate. In addition, the artificial attenuation 
mentioned earlier is applied by calling ATTENU and calls are 
made to WRITE2 or PRINT2 as required. Finally the range is 
checked to see if it has advanced to the maximum range 
specified for the problem. If it has, then IFD calls 
subroutine END which passes appropriate messages to the 
terminal and stops program execution. 

D. SUBROUTINES 

1 . Subroutine READ 

Subroutine READ is called by IFD to read input data 
from unit number NIU = 51. Input data are read in free 
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format and data are transferred back to main program IFD via 
common blocks. READ contains error checks for (1) input 
data insufficient and (2) the final depth in the sound speed 
profile unequal to the maximum depth in the water column. 
If either of these error conditions exists, READ issues an 
appropriate error message to the terminal and stops 
execution. 

2 . Subroutine SVPW 

Subroutine SVPW calculates the vertical grid spacing 
used in the water and sediment. It also calculates the 
speed of sound at each of the grid points in the water 
column. Linear interpolation is used to calculate the sound 
speed at grid points between points on the user- spec if ied 
sound speed profile. 

3. Subroutine INITAL 

Subroutine INITAL initializes constants and 
variables. If the user inputs 0.0 for the value of the 
reference sound speed then INITAL sets the reference sound 
speed c Q to the sound speed averaged over the deepest water 
column. If the user inputs 0.0 for the value of the maximum 
range step then INITAL sets the maximum range step to the 
reference wavelength, 

DR MAX = XLAMDA . 

Setting the maximum range step to the reference wavelength 
is somewhat arbitrary; however, until the actual limit on 
the range step is better understood it serves as a rough 
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rule of thumb. Finally, if the user inputs 0.0 for the 
value of the range step along a horizontal interface then 
INITAL sets the range step to half of the reference 
wavelength, 

DRLVL = 0.5 * XLAMDA. 



The default range step along a horizontal interface is half 
the default maximum range step. 

4 . Subroutine MAT CON 

Subroutine MATCON calculates constants needed to 
compute tridiagonal matrix elements. Most of the constants 
computed in MATCON have no direct physical significance but 
contribute to computational efficiency. Attenuation in both 
the water and sediment is calculated with the help of a 
complex index of refraction n, 



n = 



BETA 



(1 + i 



54.575054 



) 



or 



n 2 = 



V 


2 

1, 1 


cj 


2 BETA 


_ c j_ 


T -L 


_ c j 


27.287527 



where 



BETA = attenuation (dB/wavelength ) 
c Q = reference sound speed (m/s) 

Cj = sound speed (m/s) at point j 

54.575054 = conversion factor used in converting 

db/wavelength to nepers/meter. 
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5 . Subroutine SFIELD 



Subroutine SFIELD calculates the Gaussian starting 
field at range r = 0. This subroutine is identical with 
that of. Lee and Botseas (1982); both yield the starting 
field suggested by Brock (1978), 

U(I) = CMPLX (PR, 0.0) 

where 



PR 

ZM 

ZS 

GW 

FK 

GA 



ZM-ZS 


2 


-ZM-ZS 


GW 




GW 



depth (m) of grid point 
source depth (m) 

Gaussian width (m) (= 2/FK) 

reference wave number (1/m) (= 21Y£/c Q ) 

Gaussian amplitude [= ( 2/FK) 1 ^/GW 



6 . Subroutine WRITE1 

Subroutine WRITE 1 outputs data to a file that is used 
by the plotting routine. This output file corresponds to 
unit file number NOU = 52. 

7 . Subroutine PRINT1 

Subroutine PRINT1 outputs formatted data to a file 
that can be sent to the printer. The output file 
corresponds to unit file number NPOUT = 55. 
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8 . Subroutine NEWSEG 



Subroutine NEWSEG is called at the start of each new 
linear bottom segment. NEWSEG computes and initializes 
variables that depend on characteristics of the segment. 
One of the variables initialized is ISLOPE which is a slope 
flag having value 1 if the bottom slopes down, 2 if the 
bottom is horizontal, 3 if the bottom slopes up, 4 if the 
bottom slopes down but must be modified because the slope is 
very small, or 5 if the bottom slopes up but must be 
modified because the slope is very small. NEWSEG also 
issues error or warning messages as appropriate. 

9. Subroutine NEWMAT 

Subroutine NEWMAT calculates matrix elements for the 
X and Y matrices. The Y matrix corresponds to the range at 
which the solution field is known and the X matrix 
corresponds to the range at which the solution field is to 
be found. The Y matrix is multiplied by the known solution 
field to obtain the right-hand side column vector needed to 
solve the tridiagonal system. 

NEWMAT sets up the tridiagonal matrix system for the 
new bottom segment and than calls TRIDG to solve the system 
at the first range step. It then calls ATTENU to apply 
artificial attenuation, calls WRITE2 or PRINT2 as required, 
and finally updates the interface pointer that indicates the 
index of the grid point at the water-sediment interface. 
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1 0 . Subroutine WRITE2 



Subroutine WRITE2 is basically a continuation of 
subroutine WRITE1 . It outputs data to a file corresponding 
to unit file number NOU = 52. This is the file that is 
used by the plotting routine. At range intervals specified 
by the user, WRITE2 outputs range, depth, and the value of 
u(r, z) . 

1 1 . Subroutine PRINT2 

Subroutine PRINT2 is basically a continuation of 
subroutine PRINT1 . It outputs formatted data to a file 
corresponding to unit file number NPOUT = 55. This file can 
be sent directly to the printer. At range and depth 
intervals specified by the user PRINT2 outputs tabular 
values of transmission loss and u(r,z). 

12. Subroutine TRIDG 

Subroutine TRIDG solves a linear tridiagonal matrix 
system. TRIDG is a modified version of subroutine TRID as 
listed in Gerald (1980). The major modifications to 
subroutine TRID involved introducing arrays CTWO and CR to 
preserve the original matrix element values and to make the 
routine more efficient. Introducing the two new arrays 
requires more storage space but results in a substantial 
savings in execution time, particularly for the case of a 
horizontal interface. 
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13. Subroutine TRIDL 



Subroutine TRIDL is a modified version of subroutine 
TRIDG. TRIDL differs from TRIDG in that it does not compute 
arrays CTWO and CR but rather uses the array values 
calculated in TRIDG. 

1 4 . Subroutine DOWN 

Subroutine DOWN updates the tridiagonal matrix and 
calls subroutine RHS to update the right-hand side for the 
case of a downward sloping interface. DOWN then calls 
subroutine TRIDG to solve the tridiagonal system of 
equations and finally updates the interface pointer. 

1 5 . Subroutine UP 

Subroutine UP performs exactly the same tasks as 
subroutine DOWN, but for the case of an upward sloping 
interface. Subroutine TRIDG is again called to solve the 
system. 

1 6 . Subroutine LEVEL 

Subroutine LEVEL is called to advance the solution 
for the case of a horizontal interface. For this case the 
tridiagonal matrix elements at the advanced range need not 
be changed from the previous calculation. Therefore, LEVEL 
need only update the right-hand side by calling RHS and 
then solve the system by calling TRIDL. 

1 7 . Subroutine RHS 

Subroutine RHS computes the right-hand side of the 
tridiagonal system by multiplying tridiagonal matrix Y by 
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the known solution field U(I). The resultant right-hand 
side column vector is stored in C(I,4). 

1 8 . Subroutine SSLOPE 

Subroutine SSLOPE is called to advance the solution 
in the case where the bottom has been modified. SSLOPE 
determines which of three cases a particular section falls 
into: a level section following a level section, a level 
section following a sloping section, or a sloping section. 
For the case of a level section following a level section 
SSLOPE calls LEVEL to advance the solution. For the case of 
a level section following a sloping section SSLOPE updates 
appropriate matrix elements, calls RHS and then calls TRIDG 
to advance the solution. And in the case of a sloping 
section SSLOPE updates matrix elements and calls either DOWN 
or UP as appropriate. 

1 9 . Subroutine ATTENU 

Subroutine ATTENU applies artificial attenuation to 
the bottom portion of the sediment layer as suggested by 
Brock (1978). The artificial attenuation matrix ATT(1) is 
calculated in subroutine NEWMAT. 

20 . Subroutine END 

Subroutine END is called when the solution field has 
reached the maximum range specified. END sends appropriate 
messages to the terminal and stops execution. The messages 
are applicable to the program as installed on the NPS 
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computer but may not be appropriate for the program if 
installed on another system. 

E . INPUT DATA 

1 . Input File 

The input data must be stored in a file 
corresponding to unit number NIU as assigned in subroutine 
READ. In its present form READ sets the input unit number 
to NIU = 51. If the user prefers to read the data from a 
different unit (for example, a card reader), then variable 
NIU in READ should be set equal to the appropriate unit 
number . 

2. Input Format 

The input data is read in free format. The input 
card images (or input cards) are arranged as follows: 

CARD CONTENTS 
1 FRQ , ZS, ZR, CO, N 

where 

FRQ = frequency (Hz) 

ZS = source depth (m) 

ZR = receiver depth (m) 

(program will reset to depth of nearest 
grid point) 

CO = reference sound speed (m/s) 

If CO = 0.0, CO is set to the sound speed 
averaged over the deepest water column. 
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N = number of vertical grid points 

If the user desires that every integer 
depth value correspond to a grid point then 
(neglecting dimensions) N should be set to 
an integer multiple of ZLYR2, the depth of 
the pressure release surface. 

CARD CONTENTS 

2 RMAX , DRLVL , DRMAX , WDR , PDR , PDZ 

where 

RMAX = maximum range (m) of solution 

DRLVL = range step (m) for marching solution 
along horizontal interface 
If DRLVL = 0.0, then DRLVL is set to 1/2 
wavelength. 

If DRLVL is greater than DRMAX, then 
DRLVL is set to DRMAX. 

DRMAX = maximum allowable range step (m) 

If DRMAX = 0.0, then DRMAX is set to 1 
wavelength . 

WDR = range increment (m) at which solution is 
written to file used by plotting routine 

PDR = range increment (m) at which solution is 
printed 

PDZ = depth increment (m) rounded to nearest DZ 
at which solution is printed 
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CARD CONTENTS 



3 

4 

5 



N 

N+1 

N+2 



BR( 1 ) , BZ ( 1 ) 

BR( 2 ) , BZ ( 2 ) 

BR( 3 ) , BZ ( 3 ) 




BOTTOM PROFILE 

Range and depth of water (m). 
Maximum number of points = 100, 
Program will reset depths to 
nearest grid point. 



-1 



-1 



This card marks end of bottom 
profile. 



ZLYR1 , RH01 , BETA1 



where 



N+3 

N+4 



N+M 



ZLYR1 = maximum water depth (m) 

RH01 = density of water (g/cm^) 

BETA1 = attenuation in water (dB/meter) 

If BETA1 is less than 0.0 then program 
calculates BETA 1 with an empirical 
formula (Brock, 1978). 

ZSVP(1), CSVP(1) SOUND SPEED PROFILE 

ZSVP(2), CSVP(2) Depth (m ) and sound speed 

L (m/s ) . 

. . \ ZSVP(1) must equal 0. 

The last depth must equal 

ZLYR1 . 
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CARD CONTENTS 



N+M+1 ZLYR2, RH02, BETA2 , C2 
where 

ZLYR2 = depth (m) of pressure release surface at 
bottom of sediment layer 
RH02 = density (g/cm^) of sediment 
BETA2 = attenuation (dB/wavelength ) in sediment 
C2 = sound speed (m/s) in sediment 

N+M+2 ZABLYR 
where 

ZABLYR = depth (m) of upper surface of artificial 
attenuation layer in sediment. 

ZABLYR should be about 3/4 of ZLYR2 . 

F . PROGRAM OUTPUT 

1 . Output Printer File 

The program outputs formatted data to a file 
corresponding to unit number NPOUT which is set to 55. This 
formatted data file may be sent to the printer if desired. 
On another system the user may elect to assign NPOUT to the 
unit number corresponding to the printer and thereby send 
the formatted data directly to the printer. 

2. Output Plotter File 

The program outputs unformatted data to a file that 
is used by the plotting routine. The unit number for this 
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file is NOU which is set to 52. The output data in this 
file are stored as follows: 

LINE CONTENTS 

1 RMAX 
where 

RMAX = maximum range (m) of solution 

2 RA, ZR, U 

3 RA, ZR, U 

• • • • 

• • • • 

where 

RA = range (m) 

ZR = depth (m) 

U = complex u at specified range and depth 

3 . Terminal Output 

Certain WRITE statements in the program specify unit 
number 6. Unit number 6 on the NPS computer for an 
interactive program corresponds to terminal output. If the 
program is not run interactively then WRITE statements with 
unit number 6 may be deleted. Any pertinent information 
passed to the terminal during an interactive run is also 
passed to unit number 55. 

G. PLOTTING PROGRAM 

Appendix C contains a listing of the IFD plotting 
program installed on the NPS computer. The filename and 
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filetype of the program are PLOTIFD FORTRAN. This program 
was written separately from the IFD program because it was 
recognized that different computer installations have 
different plotting facilities. For users with different 
facilities the program will be a helpful reference. 

The program reads data from unit number NOU = 52 which 
corresponds to a file with filename/f iletype IFDOUT PLOTTER. 
Details concerning using the plotting routine are included 
in Appendix B. 
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IV. TEST CASES 



A. HORIZONTAL INTERFACE CASES 

The IFD program treats a horizontal interface using the 
same theoretical approach as the IFD program published by 
Lee and Botseas (1982). Throughout the remainder of this 
thesis the Lee and Botseas (1982) program will be called the 
FINITE-DIFFERENCE program and the program presented in this 
thesis will be called the IFD program. Two cases were run 
to confirm that the IFD program is in agreement with the 
FINITE-DIFFERENCE program for horizontal interfaces. 

1 . Isospeed Shallow Water 

This case , first published by Jensen and Kuperman 
(1979), considers propagation in a shallow water, isospeed 
environment. The water depth is 100 meters and the solution 
field is calculated out to 25 kilometers. The sound speed 
is 1500 m/s in the water and 1550 m/s in the sediment. 
Density and attenuation in the sediment are 1.2 g/cm J and 1 
dB/wavelength respectively. The source and receiver are 
both at 50 m and the source frequency is 500 Hz. 

Solutions obtained using a normal mode program 
(SNAP), a split-step fast Fourier transform program ( PAREQ) 
and the FINITE-DIFFERENCE program are shown in Figure 2. 
SNAP and PAREQ are programs that were developed at SACLANT 
Centre and are discussed further in Jensen and Kuperman 
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PROPAGATION LOSS (dB) 



SNAP 




F 

SD 

RD 



500 Hz 
50 m 
50m 






Figure 2. Propagation Loss Versus Range for Shallow Water 
Case; SNAP, PAREQ and FINITE-DIFFERENCE Results 
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(1979). The solution obtained using the IFD program is 
shown in Figure 3. All of the solutions are in excellent 
agreement . 

The input runstream that produced the results shown 
in Figure 3 for the IFD program is as follows: 

Input Runstream 
500 50 50 0 500 

25000 5 5 50 5000 50 

0 1 00 

25000 100 

-1 -1 

100 1.0 - 1.0 

0 1500 

100 1500 

250 1.2 1.0 1550 

200 



2 . Horizontal Interface 

This case, called the "horizontal interface problem" 
in Lee and Botseas (1982), considers propagation in an 
environment with the sound speed profile shown in Figure 4. 
Source frequency is 100 Hz, source depth is 30 m and receiv- 
er depth is 90 m. The density in the bottom is 2.1 g/cm^ 
and the sound speed in the bottom is 1505 m/s. No attenua- 
tion is applied in the water or sediment using complex 
indices of refraction; however, artificial attenuation is 
applied in the lower portion of the sediment layer. 
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Figure 3. Propagation Loss Versus Range for Shallow Water 
Case; IFD Results 



SOUND SPEED (m/s) 
1480 1490 1500 




Figure 4. Horizontal Interface Case 

The solution obtained using the IFD program is shown 
in Figure 5. This solution is in excellent agreement with 
the FINITE-DIFFERENCE solution shown in Lee and Botseas 
( 1 982 ) . 
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Figure 5. Propagation Loss Versus Range for Horizontal 
Interface Case; IFD Results 



The input runstream for this case is as follows: 



Input Runstream 



100 

20000 

0 

20000 

-1 

240 

0 

1 20 
240 
1 200 
51 2 



30 

2 

240 
240 
-1 
1 .0 
1500 
1 498 
1500 



90 0 

2 50 



0.0 



600 
1 0000 



2.1 0.0 1505 



50 



B. RANGE-DEPENDENT CASES 

The following range -dependent cases were solved by 
Jensen and Kuperman using SNAP, the normal mode program, and 
PAREQ, the SSFFT program (Jensen and Kuperman, 1979). The 
cases were also solved by Lee and Botseas using the FINITE- 
DIFFERENCE program (Lee and Botseas, 1982). The FINITE- 
DIFFERENCE program treats the sloping interface as a "stair 
step" and uses the interface conditions appropriate for a 
horizontal interface. The IFD program handles the sloping 
interface using the interface treatment derived by Lee and 
McDaniel (1983). 
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1 . Deep- to- Shallow Water 



This case considers propagation in an environment as 
depicted in Figure 6. The problem is solved for a bottom 
with a 8.5 degree upslope and one with a 0.85 degree 
upslope. Source frequency is 25 Hz, source depth and 
receiver depth are 25 meters. Sound speed in the water is 
1500 m/s. In the sediment, sound speed is 1600 m/s, density 
is 1.5 g/cm 3 and attenuation is 0.2 dB/wavelength . 



0 10 20 30 40 km 




The results for the 8.5 degree upslope case as 
produced by SNAP,PAREQ, and FINITE-DIFFERENCE are shown in 
Figure 7. The results as produced by IFD are shown in 
Figure 8. The difference between the results produced by 
SNAP and PAREQ is attributed to failure of the "adiabatic" 
theory underlying the SNAP program (Jensen and Kuperman, 
1979). As determined from the input runstream the FINITE- 
DIFFERENCE results were obtained using 1.0 rather than 1.5 
g/cm 3 for the density in the sediment (Lee and Botseas, 
1982). 
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Figure 7. Propagation Loss Versus Range for Deep-to-Shallow 
Water Case, 8.5 Degree Slope; SNAP, PAREQ and 
FINITE-DIFFERENCE Results 
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Figure 8. Propagation Loss Versus Range for Deep- 
to-Shallow Water Case, 8.5 degree Slope; 
IFD Results 



The input runstream for the IFD program that 
produced the results shown in Figure 8 is as follows: 

Input Runstream 



25 


25 


25 


0 


1 000 


40000 


1 0 


0 


1 00 


10000 


0 


350 








1 0000 


350 








1 2000 


50 








40000 


50 








-1 


-1 








350 


1 .0 


-1 






0 


1 500 








350 


1500 








1000 


1.5 


0.2 


1600 




750 











Appendix D contains the printed output produced by the IFD 
program using the above runstream. 

The results for the 0.85 degree upslope case as 
produced by SNAP and PAREQ are shown in Figure 9. The 
results produced by IFD are shown in Figure 10. 

2 . Shallow- to -Deep Water 

This case considers propagation in the environment 
depicted in Figure 11. The environment is exactly the same 
as the deep -to- shal lo w water environment except that the 
shallow and deep portions have been reversed and thus the 
bottom slopes down rather than up. 
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Figure 9. Propagation Loss Versus Range for Deep-to-Shallow 
Water Case, 0.85 Degree Slope; SNAP and PAREQ 
Results 



The results for the 8.5 degrees downslope case as 
produced by SNAP and PAREQ are shown in Figure 12. As 
before, the difference between SNAP and PAREQ is attributed 
to failure of the SNAP program. The results produced by IFD 
are shown in Figure 13. 

The results for the 0.85 degree downslope case are 
shown in Figures 14 and 15. 

3 . Comments 

Differences between the results obtained using the 
SNAP and PAREQ programs for the range-dependent cases are 
discussed in Jensen and Kuperman (1979). The major 
differences are attributed to the violation of the 
adiabatic assumption in the SNAP program. 
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Figure 10. Propagation Loss Versus Range for Deep- 

to-Shallow Water Case; 0.85 degree Slope; 
IFD Results 



0 




Figure 11. Shallow- to-Deep Water Case 




RANGE (km) 



Figure 12. Propagation Loss Versus Range for Shallow-to- 
Deep Water Case. 8.5 degree Slope; SNAP and 
PAREQ Results 
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Figure 13. Propagation Loss Versus Range for 

Shal low- to-Deep Water Case, 8.5 degree 
Slope; IFD Results 




Figure 14. Propagation Loss Versus Range for Shallow-to- 
Deep Water Case, 0.85 degree Slope; SNAP and 
PAREQ Results 



Figure 16 shows the results produced by IFD for the 

8.5 degree, deep-to- shal low water case if 1.0 rather than 

1.5 g/cm^ is used for the density of the sediment. The IFD 
results obtained using 1.0 g/cm^ are in very close agreement 
with the PAREQ results shown in Figure 7. However, when the 
correct value of 1.5 g/ cm J is used the slope of the 
propagation loss curve beyond 15 km is less steep (Figure 8) 
and the results do not agree as well with the results 
produced by PAREQ. The observation that the slope of the 
propagation loss curve becomes less steep when 1.5 g/cm J is 
used is qualitatively consistant because a higher density 
difference means a higher reflection coefficient which in 
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Figure 15. Propagation Loss Versus Range for 

Shallow- to-Deep Water Case, 0.85 degree 
Slope; IFD Results 
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turn results in more energy confined within the water layer. 
The difference in slope between the PAREQ and IFD results is 
attributed to the PAREQ program's apparent failure to 
account for the affects of the density discontinuity at the 
water-sediment interface. Because the IFD program correctly 
accounts for density discontinuities the results produced by 
IFD are believed to be more accurate than those of PAREQ 
when interface interaction is important. 
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V. COMMENTS AND CONCLUSIONS 



The IFD method is an efficient, stable method for 
solving the parabolic equation. Use of the IFD method is 
particularly advantageous in shallow water environments 
where the water-sediment interface is an important 
parameter. 

The IFD program presented in this thesis incorporates 
continuity of pressure and continuity of the normal 
component of particle velocity across horizontal and sloping 
interfaces. The program's capability to incorporate the 
exact interface conditions on a sloping interface, to 
automatically determine step-size, and to modify the bottom 
as required for the case of a very gently sloping bottom are 
important features. 

Projected program enhancements include wide angle 
propagation (Lee and Gilbert, 1982), range-dependent sound 
speed profiles in the water, range-dependent sound speed 
profiles in the sediment layer, and multiple sediment layers 
with horizontal or sloping interfaces. These enhancements 
are listed in their approximate order of importance. The 
program's modular construction and structured style will 
facilitate implementation of these enhancements. 
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APPENDIX A: IMPLICIT FINITE-DIFFERENCE PROGRAM LISTING 
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APPENDIX B 



RUNNING THE IMPLICIT FINITE-DIFFERENCE 
PROGRAM ON THE NPS COMPUTER 



A. INTRODUCTION 

This appendix describes one simple procedure for running 
the IFD program on the NPS computer. 



B. COPYING FILES ONTO USER'S DISK 



Four files are needed to run the IFD program; the 
filenames and filetypes are 



IFD 

IFD 

PLOTIFD 
PLOT IFD 



FORTRAN 

EXEC 

FORTRAN 

EXEC 



The files are available from a computer account maintained 
by the underwater acoustics curriculum. To link with this 
account a,nd obtain copies of the files the user should 
proceed as follows: 

(1) Log on terminal 

(2) Enter: CP LINK 01 60P 191 195 RR 

(3) When prompted for read password enter: UX 

(4) Enter: ACC 195 C 

(5) Enter: COPYIFD 

At this point the four files should reside on the user's A 
disk . 
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C. RUNNING THE IFD PROGRAM 



Before running the IFD program the user should define 
additional storage space, compile the program, and set up 
the input data file. To define additional storage space, 

(1) Enter: DEF STOR 1M 

(2) Enter: I CMS 

These two commands need only be entered one time for each 
terminal session; the storage space will remain for the 
entire session. 

To compile the program, 

Enter: FORTGI IFD 

The program need be compiled only one time unless the 
program is changed, in which case the new version should be 
recompiled . 

The final step before running the program is setting up 
the input data file which has filename and filetype IFDIN 
DATAIN. The user must create or modify this file so that it 
contains input data as described in Section III.E.2 of this 
thesis. For more information concerning how to create or 
modify files see NPS Technical Note TN-VM-02 which is 
available in the computer consultant's office. 

If the above steps are accomplished, the user can then 
run the program with 

Enter: IFD 

Shortly after entering this command the user will be 
prompted for a run identification. The run identification 



is an arbitrary identification label that will appear on the 
output printer file. Enter run identification as desired. 

At the end of the run the user is informed of the two 
output data files generated by the program. If desired by 
the user the output printer file (IFDOUT PRINTER) may be 
sent to the printer. The output plotter file, IFDOUT 
PLOTTER, serves as an input file for the plotting program. 

D. RUNNING THE PLOTTING PROGRAM 

A Tektronix-61 8 terminal is used to run the plotting 
program. The first step is to log onto the terminal in the 
normal manner and then define additional storage space by, 

(1) Enter: DEF STOR 1M 

(2) Enter: I CMS 

The plotting program has filename and filetype PLOTIFD 
FORTRAN. To compile the program. 

Enter: FORTGI PLOTIFD 

Unless the program is changed it need only be compiled one 
time. To run the program 

Enter: PLOTIFD 

The user will be prompted for axes and smoothing 
information; enter values and responses as appropriate. The 
transmission loss curve will be displayed on the CRT screen, 
a hard copy may be obtained by pushing the HARD COPY button 
under the screen, and the screen may be cleared by pushing 
the ENTER key. 
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APPENDIX C: PLOTIFD PROGRAM LISTING 



OOOOOOOOOOOOOOOOOOOOO 3000000000000000 OOOOOOOOOOO 

-4r\j o 't m so rs. oo 0>0 ^ vf so ao on o .-4 <M (o vfio so ^ ^0 O' O rg O ^ oo o — • (\j o co 

O O o O O O o o Or-4 rH r-H r-4 -H ^ r-1 r-4 rH (\4 C\| og (M (NJ CM CM c\j i\j rg rn rn rn O O rn O ro O CO 'T nT -r *>4“ -r 'T 
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APPENDIX D 



EXAMPLE OF IMPLICIT FINITE -DIFFERENCE PRINTED OUTPUT 
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